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Meteor science 

Estimating meteor rates using Bayesian inference 

Geert BarentseT^ Rainer Arlt^ Hans-Erich Frdhlich^ 

A method for estimating the true meteor rate A from a small number of observed meteors n is derived. We 
employ Bayesian inference with a Poissonian HkeHhood function. We discuss the choice of a suitable prior and 
propose the adoption of Jeffreys prior, P{X) = A~°'^, which yields an expectation value E{\) = n + 0.5 for any 
n > 0. We update the ZHR meteor activity formula accordingly, and explain how 68%- and 95%-confidence 
intervals can be computed. 



1 Introduction 

The formula commonly used to estimate the Zenithal 
Hourly Rate (ZHR) of meteors as given in the Hand- 
book of the International Meteor Organization (IMO; 
Rendtel & Arlt, 2008) and used in many meteor activ- 
ity graphs (e.g., Arlt & Barentsen, 2006), is given by: 



E(ZHR) = "^"L"- ' with a = + ' 



T 



T 



(1) 



where ntot is the total number of meteors counted in a 
number of observing intervals (ntot = X^i ^i) ^"^^ 
the observing time weighted by a given correction factor 
for each interval (T = Tcs^i/Ci). 

The use of ntot + 1 rather than ntot often surprises 
observers, because it yields a ZHR which is larger than 
zero even when no meteors are observed. Whilst Arlt 
(1999) already indicated that ritot + l is used to account 
for the effects of small-number statistics, this paper will 
explain the formula in more detail. However, contrary 
to the earlier publications, we will suggest that ntot+0.5 
rather than ntot -l- 1 is the most appropriate formula. 

We will first describe the problem of small-number 
statistics in §2. We then describe the solution using 
Bayesian inference in §3, followed by a discussion on the 
choice of the prior assumptions in §4. In §5 we explain 
how confidence intervals may be computed, and in §6 we 
provide examples. Finally in §7 we discuss the effect of 
correction factors and in §8 we present the conclusions. 

2 Problem description 

As explained by Bias (2011) in a recent issue of this 
journal, the observed meteor rate often tends to un- 
derestimate the true meteor rate due to small-number 
statistics. This may be understood using the follow- 
ing example: if we would attempt to estimate the fre- 
quency of winning the lottery based on a small group of 
10 players, we are most likely to find that none of these 
players have ever won and that the winning frequency 
equals zero. Of course the true probability of winning 
the lottery is slightly larger than zero, but the quan- 
tity cannot reliably be obtained by extrapolating from 
a small number of players. 

An identical effect occurs when the ZHR is extrapo- 
lated from a low number of meteors. Indeed, even when 



zero meteors are observed in an interval of finite length, 
the true average rate may well have been larger than 
zero because we know that the interplanetary space is 
not empty. The observer might have been unlucky, or 
the interval might have been short with respect to the 
true rate. Such situation may occur even during major 
showers, e.g. when computing rates for 1-minute inter- 
vals. The ZHR formula must take this into account in 
order to produce reliable estimates in all situations. 

3 Solution: Bayesian inference 

A common technique used in statistics to estimate pa- 
rameters in the presence of sparse data is called Bayesian 
inference. In brief, one constructs a parameterized model 
which one thinks describes the source of the data. The 
probability of this model to have produced the data 
is then computed for each possible set of parameters, 
taking into account any known prior constraints on the 
parameters. The resulting set of probabilities is called 
the posterior distribution, from which expectation val- 
ues and confidence intervals for the free parameters may 
be derived. 

In our case, the data is the observed meteor rate 
n (in arbitrary time T). Our only free model param- 
eter is the true meteor rate A (i.e. ZHR). There are 
many values of A which may explain a given n, each 
having the conditional probability P{X\n). This nota- 
tion means the probability of finding A given that one 
has seen n meteors. This is the posterior probability 
described earlier which may be estimated using the the- 
orem by Bayes: 



P{X\n) 



P{n\\)P{\) 
1^) ' 



(2) 



where P{n\\) is the generative model, i.e. the known 
probability to see n meteors for a given true rate A. 
We may assume that meteors appear in a random way 
following a Poissonian law for independent events: 



P{n\\) 



n ^— A 



A"e 



(3) 



The function P{n) serves to normalize the distribution 
to unity, while P(A) expresses what we know about how 
probable each of the possible true rates A are. This 
function is called the prior and, ideally, should be a 
probability distribution on its own. The challenge is to 
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decide what we can assume about P(A) beforehand? 



= n+ 1. 



(5) 



4 Choice of the prior 

The criteria for choosing a suitable prior P(A) is a sub- 
ject of debate in the statistical community. On one 
hand, one may decide to construct a prior based on pre- 
vious evidence, for example the meteor activity in the 
past decade. This is called an informative prior. On the 
other hand, one may prefer a prior which contains only 
vague or general information and is not biased towards 
past observations. This is called an objective prior. 

Whether or not it is appropriate to include previous 
observations in the computation of meteor rates is a 
philosophical question. However, given the intrinsically 
variable nature of meteor showers, we suggest that an 
objective prior is the only practical approach. We dis- 
cuss the choice of such prior in what follows. 

4.1 Uniform prior 

A natural choice is a uniform prior P{X) = constant, 
i.e. all values of A are assumed to be equally likely. A 
uniform prior leads to a distribution which cannot be 
normalized, and is called an improper prior. Any limit 
on ZHR, and be it very high, makes the distribution 
normalizable though, that is, asymptotically a uniform 
prior is not a problem. 

Let us see what the uniform prior implies for the 
inferred rate. The expectation value of the true meteor 
rate is defined as the sum of every possible value of A 
times its posterior probability. (In this paper, we also 
use the terms "average" or "mean" as synonyms for 
expectation value.) For a continuous quantity such as 
A, the expectation value is a normalized integral over 
all A which, for a uniform prior, reads: 



The expectation value of the true activity rate is 
thus given by (n + 1) with spread of the distribution 
of (7 = \/n + 1. These quantities arc then simply di- 
vided by the weighted time T to obtain the ZHR. This 
explains the formula given in the IMO Handbook. 

4.2 Exponential prior 

A uniform prior may not be ideal if we want to express 
the fact that very large rates are very unlikely. Almost 
all rates ever obtained are below say 1000 meteors per 
hour. A suitable prior will decrease with rate, and it is 
mathematically convenient to use power functions like 



P(A) = 1/Ai 



(6) 



where < a < 1. Such a power-law has the advantage 
that we can use F-functions for the derivation of -E(A). 
The resulting expectation value now involves the prior 
and is 



J XP{X)P{n\X)dX 

E{X) = ^ , 

J P{X)P{n\X)dX 





(7) 



where the lower integral is to normalize the posterior 
distribution. Now, let's make use of the F-functions for 
which 



n! = nF(n) = T{n + l) = j X^-e'^dX 





(8) 



holds. We do not need to know how F is computed, we 
only need its properties to derive the expectation value: 



E(A) = 



AP(A|n) dA 
A^e^^dA 



f° 
Jo 

n + 1. 



A"+i 
(n+1)! 



-^dA 



(4) 



Similarly, the error estimate for E(A) follows from the 
mathematical definition of the variance, often referred 
to as a^: 



a^X) = 



E[(A-E(A))2] 

/>00 

/ P(A|n)(A-(n-|-l))2dA 
Jo 



-I 



(n+ l)(n + 2)A"+2e-^ 
(n + 2)! 



+ 



2(n + l)^A"+ie-^ 

{n + l)\ 
{n+lfX^e-^ 



+ 



dX 



(n+l)(n + 2)-(n+l)^ 



A" A ^-A 



E{X) 



« r(n-l-l) 



e-^dX 





r(n+l+a) 

r(n+i) 

r(n+a) 

{n + a)T{n + a) 

r(n + a) 
n + a. 



Similarly, the variance is: 



(9) 



(t2(A) = E (^X-E{X)y 

/P(A)P(n|A)(A- (n + a)) dX 

_ _^ ' 

oo 

J P{X)P{n\X)dX 



_ {n + l + a){n + a)r{n + a) 
~ F(n + a) 

2(n + a)^r(n + a) + {n + a)^F(n + a) 
F(n + a) 
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= {n+l + a){n + a)-2{n + af + {n + af 
= n + a, (10) 

The expectation value of the true activity rate is 

thus given by (n + a.) with a spread cr = ^/n + a, for 
a given prior 1/A^~". Again, these priors are normahz- 
able for < a < 1 by enforcing a limitation of the vahd 
range. 

One question remains: which is the most appropri- 
ate value to adopt for a? On one hand, the uniform 
prior (a = 1) yields the expectation vahie n + 1, which 
is Ukely to overestimate the meteor rate due to assump- 
tion that any arbitrarily large rate is equally likely as 
the zero rate. On the other hand, the prior 1/A (a = 0) 
yields a "traditional" extrapolation E{X) = n, which is 
likely to underestimate the rate as explained previously. 
It appears appropriate to adopt a prior somewhere in- 
between < Q! < 1. 

4.3 Jeffreys prior; a = 0.5 

The problem of choosing a suitable prior for a Pois- 
sonian process exists in other fields (e.g. radioactive 
decay counts, neutrino detections). An axiomatic solu- 
tion has previously been proposed by Harold Jeffreys. 
He required that a prior should be "invariant under 
reparameterization" , i.e. a prior should not depend on 
the variable investigated (Jeffreys 1946, 1961; Kass & 
Wasserman 1996). One could, for example, be inter- 
ested in the rate A as well as in the mean time between 
events /z = 1/A. The priors for both expectation values 
should be compatible. For general relations between 
different quantities, the requirement of compatibility is 
written mathematically as 



P{\)d\ = P{ii)dn. 



(11) 



For a Poissonian distribution, such a general compati- 
bility is achieved for a prior a = 0.5. In this case, the 
estimate is not specific to either computing the rate or 
the mean time lapse or something else. 

Because there are no further statistical principles to 
decide which prior is to be preferred, we suggest the use 
of Jeffreys prior 



P(A) = 1/A' 



0.5 



(12) 



for future estimates of the rate. The expectation value 
for the meteor rate is then 



£;(A) = (n + 0.5) ± Vn + 0.5. 



(13) 



An illustration of the posterior distribution P{\\n) un- 
der the assumption of Jeffreys' prior is shown in Figure 1 
for different values of n. 

Note that in most cases, when n is large (n ^ 0.5), 
the differences between the various priors are negligable. 

5 Confidence intervals 

The standard deviation a = ^/n-\- 0.5 characterizes the 
spread in the posterior distribution around the expec- 
tation value. In the case of a Gaussian distribution, 



the standard deviation corresponds to a confidence in- 
terval (i.e. the 68%-confidcncc interval of a Gaussian 
distribution is located between —la and +lo" from the 
mean). 

However, the posterior P{\\n) docs not follow a 
Gaussian shape and is more similar to a Poissonian 
distribution (though a Gaussian shape is approached 
for large n). The posterior is asymmetric with a tail 
towards high meteor rates, which makes it somewhat 
misleading to characterize the uncertainty with a single 
number. A better way to characterize the uncertainty 
is to compute the (asymmetric) error margins of the 
68%-confidence interval. This may be done as follows. 

Given the posterior distribution 



P{\\n) 



(14) 



T{n + a) 

The corresponding cumulative distribution function is 



-P'(A < Amarg|n) = 



-^marg ^n—l + a^ — t 



r{n + a) 



dt 



^ Tia + n) ' ^^^> 



with the incomplete F function which for integers n > 
can be computed as 



r(n,A^ 



marg/ 



(n- l)!e- 



n-l ^fc 



k=0 



k\ 



(16) 



By integrating the cumulative distribution function nu- 
merically for different probabilities (P' = 2.5%, 16%, 

84%, and 97.5%), we obtain useful quantiles which cor- 
respond to the central 68%- and 95%-confidence inter- 
vals. 

These quantiles are shown in Table 1, given as a mul- 
tiplier of the true rate. For example, after computing 
the ZHR, one may look up the corresponding relative 
margins for a given ntot in Table 1 and obtain the 68%- 
interval by computing ZHR • Jes.iow (negative margin) 
and ZHR - ^es.high (positive margin). 

It is interesting to note that the 68% interval ap- 
proaches symmetry from ntot ^ 3, while the 95% inter- 
val is more sensitive to the wings and remains asym- 
metric even beyond ntot ^ 1000. 

For ritot larger than about 30, the 68%- interval ap- 
proaches a Gaussian shape and the margins can be ap- 
proximated using (T = v^ntot + 0.5/r (= ZHR/\/ntot + 0.5). 

Finally, we remind the reader that the margins in 
Table 1 only represent the uncertainty which is due to 
Poissonian statistics. The true uncertainty of a ZHR es- 
timate is likely to be somewhat larger due to observing 
errors (e.g., uncertainty in the limiting magnitude de- 
termination) . Fortunately, the impact of such observing 
errors is likely to be small when data from a sufficient 
number of independent observers is averaged. 

6 Examples 

Now let us consider the formula using a few examples. 
If one meteor was seen in five minutes, the rate equals 
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Figure 1: Posterior distribution P{X\n) of the true meteor rate under Jeffreys' prior (a = 0.5), plotted for different 
observed rates (n = 0, 3, 6). Vertical lines indicate the position of the expectation values E(A) = n + 0.5. 



Table 1: Margins of the 68%- and 95%-confidence in- 
tervals, given as multipliers to the ZHR. After comput- 
ing ZHR = (ritot + 0.5)/r, the zb-values can be ob- 
tained by computing ZHR • Jes.iow (negative margin) 
and ZHR • (568, high (positive margin) for the appropriate 
ntot- 



ntot 



395, low 



068. low 



J68,high 



395, high 






-1.00 


-0.96 +0.99 


+4.02 


1 


-0.93 


-0.72 +0.73 


+2.12 


2 


-0.83 


-0.59 +0.59 


+1.57 


3 


-0.76 


±0.51 


+1.29 


4 


-0.70 


±0.45 


+1.11 


5 


-0.65 


±0.41 


+0.99 


6 


-0.61 


±0.38 


+0.90 


7 


-0.58 


±0.36 


+0.83 


8 


-0.56 


±0.34 


+0.78 


9 


-0.53 


±0.32 


+0.73 


10 


-0.51 


±0.30 


+0.69 


11 


-0.49 


±0.29 


+0.66 


12 


-0.48 


±0.28 


+0.63 


13 


-0.46 


±0.27 


+0.60 


14 


-0.45 


±0.26 


+0.58 


15 


-0.43 


±0.25 


+0.56 


16 


-0.42 


±0.24 


+0.54 


17 


-0.41 


±0.24 


+0.52 


18 


-0.40 


±0.23 


+0.50 


19 


-0.39 


±0.22 


+0.49 


20 


-0.39 


±0.22 


+0.48 


22 


-0.37 


±0.21 


+0.45 


24 


-0.36 


±0.20 


+0.43 


26 


-0.34 


±0.19 


+0.42 


28 


-0.33 


±0.19 


+0.40 


30 


-0.32 


±0.18 


+0.38 



E(ZHR) = 



i8.ol}^:J. 



When zero meteors are 
seen during five minutes, the rate equals ZHR — 6.0^5 g. 
The error margins are actually so close to symmetric 
that we can always give a single value for the error bars. 



i.e. ZHR = 18+13 and ZHR = 6 ± 6 respectively 
Although the rounded margins suggest so, the lower 
margin is not zero! 

If zero meteors were seen in four hours, the rate 
equals to ZHR = 0.125^q |^2o Indeed, rates can only be 
constrained to values close to zero when no meteors are 
observed for a very long period. A more "normal" case 
for a minor shower would be say a total of 12 meteors 
in a total of 4 hours, delivering ZHR = 3.1 ± 0.9. A 
larger number of meteors of say 34 meteors in 11 hours 
gives simply ZHR = 3.1 ± 0.5 with the above error for 
large meteor numbers. 

7 The influence of correction factors 

In the previous sections we have ignored the specific 
correction factors which are used to obtain a standard- 
ized ZHR (e.g. to account for limiting magnitude and 
radiant elevation). Given two rates; one without cor- 
rection, A and one with correction. A, we have assumed 
there is a factor / which does not depend on the rate: 



/A- A 



(17) 



Indeed, for a set of N observing periods being combined 
in one expectation value for A, all having different cor- 
rection factors fi, we obtain 



/ 1 A (/^^)'"(/?^r-({"^)"" e-i:/.^dA 

J A-'- " ni!n2!---njv! 



f 1 J/i^)''^(^A)"2...(/„A)--» ^,^^,,^^ 



r(ntot + l + a) (E/.)""^^" 

(E /»)"'°'^'^" r(ntot + «) 

»tot + a 



(18) 



In other words, our method may be applied regardless 
of the values of the correction factors. 
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8 Discussion and conclusion [8] Rendtel, J., Arlt, R. (eds.) (2008) "IMO Handbook 

In conclusion, we recommend to compute ZHR values ^or Meteor Observers". IMO, Potsdam, 

using the term n + 0.5: 

E(ZHR). (19) 

Jeff sm /iR 

with error margins: 

AZHR=^^. (20) 
V"tot + 0.5 

Note that for ntot < 30, the error margins listed in 
Table 1 should be used instead of Eqn. [20] to obtain a 
68%-confidence interval. 

This method of computing the ZHR adds a small 
bias taking into account the asymmetry of possible rates. 
In particular, it is essential to adopt the method when- 
ever rates are based on less than ~ 10 meteors. Such sit- 
uations commonly occur when a major shower is anal- 
ysed using very short (e.g. 1-minute) intervals. When 
computing a rate based on a number of observing peri- 
ods (indexed with i), never ever compute the rate from 
rii + 0.5 for individual periods and average them after- 
wards. A more accurate estimate is based on the sum 
of meteors from these observing periods, and hence on 
«tot +0.5. Finally, it is, of course, much better to have 
enough observations and large enough ritot that the sub- 
tleties of choosing a good prior are no longer important, 
i.e. ntot > 0.5. 

Comprehensive information about Bayesian infer- 
ence in general and the choice of priors can be found 
in e.g. Kass feWasserman (1996) and Bolstad (2007; 
Chapter 10 for priors). 
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